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Computational  Investigation  of  Combustion  Instabilities  in 
a  Laboratory-Scale  LDI  Gas  Turbine  Engine 
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Maurice  J.  Zucrow  Laboratories,  Purdue  University,  West  Lafayette,  IN,  47906 

and 

Venkateswaran  Sankaran5 

Air  Force  Research  Laboratory,  Edwards  AFB,  CA,  93524 


Self-excited  combustion  instabilities  in  a  lean  direct  injection  (LDI)  gas  turbine 
combustor  are  computationally  investigated.  The  model  LDI  combustor  under  study  was 
developed  to  produce  combustion  dynamics  on  demand  using  a  single  LDI  element  in  an 
axisymmetric  combustor.  Three  simulation  cases  for  this  combustor  are  considered:  non¬ 
reacting  and  reacting  flow  in  an  acoustically  open  combustor  and  reacting  flow  in  an 
acoustically  closed  combustor.  We  studied  the  dynamic  flow  features  in  the  LDI  combustor 
for  both  cases  and  investigated  the  important  modes  using  a  dynamic  mode  decomposition 
method.  Precessing  vortex  core  (PVC)  instabilities  are  indicated  as  the  critical 
hydrodynamic  mode  and  lead  to  strong  spray  and  flame  response.  In  the  acoustically  close 
chamber  simulation,  we  were  able  to  capture  self-excited  combustion  instabilities  and  the 
dominant  modes  from  simulations,  which  qualitatively  agree  with  the  experimental  results. 
The  appearance  of  pressure  peaks  in  both  simulation  and  experiment  at  about  1400  Hz 
corresponding  to  the  4L  acoustic  mode  and  at  6000  Hz  are  explained  by  the  nonlinear 
coupling  of  the  PVC  and  acoustics  modes  and  the  associated  feedback  loop. 


I.  Introduction 

Lean  Direct  Injection  (LDI)  is  a  promising  design  concept  in  gas  turbine  engines.  It  was  developed  as  a  low  NOx 
alternative  to  Lean  Prevaporized  Premixed  (LPP)  combustion  for  aircraft  gas  turbines.  In  the  LPP  design,  the 
fuel  droplets  are  atomized  and  pre-mixed  with  the  air  in  a  tube  upstream  of  the  combustor.  Although  this  design 
can  provide  low  NOx  emissions  due  to  the  uniformly  low  equivalence  ratio,  the  inherent  auto-ignition  and  flashback 
are  risks  that  may  be  too  high  for  flight  applications.  On  the  other  hand,  in  LDI  combustors,  the  liquid  fuel  is  rapidly 
dispersed  by  a  highly  inertial  air  flow  and  directly  injected  into  the  flame  zone.  While  the  issues  of  concern  in  the 
LPP  design  are  nearly  eliminated,  NOx  performance  is  compromised  due  to  the  absence  of  the  premixing  and 
prevaporization  processes.  The  key  for  improved  NOx  performance  in  a  LDI  design  therefore  lies  in  good 
atomization  and  quick  and  uniform  mixing. 

Combustion  dynamics  introduce  additional  technical  challenges  to  the  design  of  low  NOx  combustors,  and  the 
vulnerability  of  the  LDI  design  to  self-excited  pressure  oscillations  is  well  documented.  [1,  2]  As  the  combustion 
occurs  under  a  lean  premixed  condition,  the  flame  surface  can  be  more  responsive  to  acoustic  fluctuations  than  a 
typical  diffusion  flame.  The  sensitivity  of  the  chemical  kinetics  is  also  influenced  by  operation  at  low  equivalence 
ratios.  Lastly,  the  acoustic  waves  may  be  less  damped  due  to  the  absence  of  significant  vorticity  generation  across 
the  orifices  in  the  LDI  design. 
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This  study  is  aimed  at  analyzing  combustion  dynamics  in  a  laboratory  LDI  combustor  and  enhancing  our 
understanding  of  the  underlying  driving  mechanisms  by  means  of  a  high-fidelity  simulations  and  companion 
experimental  tests.  A  single-element  LDI  model  combustor  which  generates  self-excited  pressure  oscillations  was 
developed  for  the  present  study.  The  computational  simulation  considers  two  different  geometries:  acoustically-open 
and  acoustically-closed  combustors.  The  acoustically-open  combustor  allows  the  study  of  inherent  hydrodynamics 
and  heat  release  without  the  influence  of  acoustics.  In  addition  to  reacting  flow  simulations  in  these  two  combustors, 
we  non-reacting  flow  in  the  acoustically-open  geometry  is  also  studied  in  order  to  understand  the  change  of  flow 
characteristics  due  to  chemical  reactions  and  longitudinal  acoustics.  The  major  mode  shapes  due  to  the 
hydrodynamic  instabilities  are  investigated  using  Dynamic  Mode  Decomposition  (DMD).  Lastly,  the  combustion 
dynamics  in  the  LDI  combustor  are  estimated  by  the  computational  simulation  and  a  physical  interpretation  of  the 
simulation  results  is  provided  to  help  understand  the  underlying  combustion  dynamics  mechanisms. 


II.  Computational  Model  Description 


A.  Computational  Framework  (GEMS) 

The  computational  platform  for  the  present  simulations  is  the  in-house  research  code  GEMS  (General  Equation 
and  Mesh  Solver). [3,  4]  GEMS  is  a  fully  unstructured,  density-based  finite  volme  solver  with  a  second-order 
numerical  scheme  and  an  implicit,  dual  time  procedure  for  time -accuracy.  The  capabilities  of  the  code  for  capturing 
the  combustion  dynamics  and  predicting  self-excited  instabilities  have  been  successfully  demonstrated  for  rocket 
engine  combustors.  [5-7]  GEMS  solves  the  Navier-Stokes  equations  in  Detached  Eddy  Simulations  (DES)  mode 
along  with  the  continuity,  energy  and  species  equation  described  below. 
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The  quantities;  p,  V  and  p  represent  the  density,  velocity  vector  and  pressure  respectively;  h°  is  the  stagnation 
enthalpy  and  Y  is  a  species  mass  fraction  vector;  K  represents  the  turbulence  variable  vector  which  includes 
turbulence  kinetic  energy  and  specific  dissipation,  k  and  w;  In  the  viscous  flux,  D  is  the  molecular  diffusion 
coefficient;  r  is  the  stress  tensor  and  q  is  the  heat  flux;  w  and  sK  are  the  reaction  rate  and  sources  of  the  turbulence 
transport  equations  respectively.  A  pseudo-time  term  expressed  in  terms  of  the  primitive  variables, 

Qp  =  [p  Y  V  T  K]t  and  a  preconditioning  matrix,  r  ,  is  added  to  Eq.  (1),  so  that  the  equation  becomes: 


dr 
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The  matrix,  r  ,  is  chosen  to  control  the  artificial  disspation  in  the  spatial  discretization  and  the  convergence  of  the 
pseudo-time  iterations.  The  preconditioning  matrix,  T  ,  in  the  pseudo  time  term  in  Eq.  (3)  is  defined  by  starting 
from  the  Jacobian  of  the  conservative  variables  with  respect  to  the  primitive  variables,  dQ  /  dQ  ,  as  shown  below: 
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The  flux  formulation  for  a  generalized  upwind  finite  volume  approach  can  be  interpreted  as  the  average  of  the 
fluxes  on  either  side  of  the  cell  interfaces  augmented  by  an  artificial  dissipation  term, 
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F=l-(FL+FR)-l-\A\(QR-QL)=X-(FL+FR)  -  \r~'Ap \(QpR-  QpL  )  (5) 

where  the  subscripts  R  and  L  represent  values  on  the  right  and  left  side  of  the  cell  face.  The  numerical  procedure 
uses  a  second-order  approximate  Riemann  solver  to  evaluate  the  spatial  fluxes  at  cell  faces.  Second-order  temporal 
accuracy  is  achieved  by  means  of  an  implicit  dual  time  procedure  that  eliminates  factorization  and  linearization 
errors. 


B.  Turbulence  Model 

The  present  simulations  describe  the  large-scale,  time-dependent  turbulence  motion  by  means  of  a  Detached  Eddy 
Simulation  (DES)  model. [8,  9]  DES  is  a  hybrid  RANS/LES  approach  which  combines  the  RANS  method  in  the 
attached  boundary  layers  with  the  LES  method  for  the  large  separated  and  wake  regions.  A  DES  model  can  be 
obtained  from  any  RANS  model  based  on  a  local  length-scale  definition  and  our  DES  formulation  uses  the  k-omega 
two  equation  model  with  the  appropriate  modification.  This  modification  relies  on  the  local  turbulence  length  scale 
and  the  local  grid  dimensions  to  enable  the  switching  from  the  RANS  to  LES  mode.  It  activates  the  RANS  mode 
near  the  solid  surfaces  where  the  flow  is  attached,  whereas  it  invokes  the  LES  mode  in  separated  flow  regions  by 
reducing  the  dissipation  term  in  turbulent  kinetic  energy  transport  equations.  Specifically,  the  DES  model  replaces 
the  length  scale  with  the  minimum  of  the  length  scale  defined  in  turbulence  model  for  RANS  and  maximum  spacing 
of  the  local  grid  as: 

l DES  =  m^n  ( ' h-a>  ’  CDES  &  )  (6) 

where  8  =  max  (Ax,  Ay,  Az)  represents  the  maximum  grid  spacing  in  any  direction  and  CDES  is  a  model  constant  and 


set  equal  to  0.78  as  recommended  by  Travin  et  al.;  lk_(jJ  is  the  length  scale  of  Wilcox’s  k-w  two  equation  model  for 
the  turbulence  closure  in  RANS  and  it  is  defined  as 


^ ' k-co 
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where  P*  is  a  model  constant.  This  definition  of  the  length  scale  by  Eq.  (6)  ensures  the  RANS  mode  is  utilized  near 


the  wall  surface  in  which  the  high  grid  apect  ratio  is  typically  expected.  Finally,  the  dissipation  term  in  the  transport 
equation  of  the  turbulence  kinetic  energy  is  replaced  by 

okm 

P*pkw  =  ^~  (8) 

DES 

This  modification  ensures  that  the  resulting  sub-grid  model  reduces  to  a  Smagorinky-like  model  at  equilibrium. 


C.  Combustion  Model 

The  two-step  global  reduced  mechanism  by  Westbrook  is  considered  in  the  present  study.  [10]  All  six  chemical 
species,  Ci2H23,  02,  C02,  H20,  CO  and  N2,  are  solved  directly  and  turbulence  transport  in  the  species  equations  is 
modeled  using  the  classical  gradient  model  with  a  constant  Schmidt  number.  Flamelet  and  transported  FMDF 
models  are  popular  in  unsteady  LES  simulations,  but  their  capabilities  for  combustion/acoustic  interaction  problems 
are  not  well-established.  Laminar  finite  rate  chemistry  model  is  therefore  used  here,  which  solves  transport 
equations  for  the  chemical  species,  and  uses  Arrhenius  rate  expressions  to  describe  the  reaction  rate.  No  additional 
turbulence-chemistry  interactions  are  considered. 

We  consider  the  fuel  as  kerosene  (Ci2H23)  to  approximate  the  experimental  Jet- A  fuel.  A  simplified  two-step,  five 
species  global  reduced  mechanism  has  been  incorporated  in  the  present  study. 

2C12H23  +3502 - >24C0  +  23H20 

2CO  +  02  < - >2C02 

The  kinetics  of  the  two-step  global  reactions  are  represented  by  an  Arrhenius  function  of  the  form  as: 

co  =  AT"  (-Ea  / i?„r)[XA  Y  [XB  f  (9) 

where  w  is  the  production  rate,  A  is  the  pre-exponential  constant,  Ea  is  the  activation  energy,  and  n ,  a  and  b  are 
exponents,  respectively.  For  the  first  Ci2H23  oxidation  step,  the  production  rate  is  expressed  with  constants  as 
2.643xl09exp(-1.5108xl04/T)[Ci2H23]°'25[O2]L5.  For  the  second  step,  it  is  expressed  as  2.23 87x10 12exp(- 
2.0 1 43 x  1 04/T) [CO] [02]° 25 [H20]°'5  for  the  forward  CO  oxidation  reaction,  and  as  5.0xl08exp(-2.0143xl04/T)[C02] 
for  the  reverse  CO  dissociation  reaction. 
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D.  Lagrangian  Phase  Equations 

The  Lagrangian  formulation  for  the  liquid  particle  motion  is  expressed  as  a  set  of  ordinary  differential  equations 
for  the  Larangian  solution  variables,  QL ,  as: 


dQi 

dt 
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dt 
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m i 
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where  X,  U,  Th  mi  and  r  are  the  components  of  the  Largrangian  solution  vector  and  these  variables  represent  the 
position  and  velocity  vector,  temperature,  mass  and  radius  of  the  liquid  particle  respectively.  The  first  and  second 
rows  represent  the  motion  of  particle  droplets,  where  the  drag  is  considered  as  an  external  force.  In  the  source  of 


3  c  p 

second  row,  FD  is  a  drag  parameter  defined  as  FD  = - — — |U  -  V|  in  which  V  is  the  carrier  gas  velocity  vector  and 

8  r  pi 

the  drag  coefficient  follows  Putnam’s  formulation  based  on  a  spherical  drop  assumption.  [11] 
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Here,  the  Reynolds  number,  Re,  is  defined  as 


U-  V|  .  The  third  row  in  Eq.  (10)  resprents  the  energy  equation 


for  the  liquid  particle.  This  equation  assumes  a  uniform  temperature  inside  the  liquid  particle.  The  source  term 
consists  of  the  net  sensible  enthalpy  transfer,  cp(Tg  -7]) ,  and  the  latent  heat  of  vaporization,  Lv ,  and  mlCl 

represents  the  heat  capacity  of  the  liquid  particle.  The  fourth  row  accounts  for  the  mass  balance  during  vaporization. 
The  corresponding  source  term  is  the  vaporization  rate,  mv ,  which  is  evaluated  by  means  of  the  classical  D2-law: 


mv  =-2?rrpgDsSh\n(l  +  BM)  (12) 

where  Ds  is  a  diffusion  coefficient  at  the  liquid  surface  and  Sh  is  a  Sherwood  number;  BM  is  the  Spalding  mass 
transfer  number  [12]  given  by 


bm=- 


'~YF,s 


(13) 


where  YF  is  the  fuel  species,  subscripts  s  and  oo  indicate  the  surface  and  far  field  respectively.  The  last  row  accounts 
for  the  liquid  particle  radius  reduction  due  to  vaporization. 

The  set  of  Lagrangian  phase  equations,  Eq.  (10),  is  integrated  explicitly  and  the  Lagrangian  time  step  is 
determined  by  consideration  of  the  time  scales  of  the  particle  lifetime,  atomization  initiation  time  and  Eulerian  time 
step,  etc.  All  other  physical  properties  of  the  liquid  particle  are  updated  every  time  step  as  well. 


E.  Injection  and  Atomization  Model 

The  atomizer  used  in  this  study  is  a  pressure  swirl  atomizer.  The  hollow  cone  fuel  spray  is  described  by  the 
injection  of  a  series  of  discrete  liquid  parcels  which  contain  a  bunch  of  drops.  This  parcel  (or  “blob”)  is  the  particle 
solved  by  the  Lagrangian  phase  equations.  Note  the  parcel  is  not  an  actual  drop,  but  rather  a  cloud  which  contains  a 
number  of  drops.  The  parcel  undergoes  a  series  of  complicated  physical  processes,  which  are  described  by 
corresponding  submodels  as  shown  in  Table  1.  As  shown  in  Figure  1,  three  main  events  are  considered  in  the 
present  model:  free  surface  flow,  primary  atomization,  and  secondary  atomization.  Before  forming  the  parcel,  there 
is  a  free  surface  flow  inside  of  the  pressure  swirl  atomizer.  Depending  on  the  flow  conditions  and  atomizer 
geometry,  the  liquid  sheet  thickness  at  the  atomizer  exit  can  be  determined  by  Ibrahim  and  Jog’s  numerical  model 
[13].  The  liquid  sheet  around  the  atomizer  exit  disintegrates  into  ligaments  and  is  subsequently  broken  into  multiple 
droplets.  This  primary  atomization  process  is  described  by  Senecal  et  al’s  Linearized  Instability  Sheet  Instability 
(LISA)  model  [14].  Furthermore,  droplets  from  the  primary  atomization  encounter  aerodynamic  forces  due  to  the 


4 

American  Institute  of  Aeronautics  and  Astronautics 


strong  swirling  air  flow,  leading  to  variety  of  secondary  atomization  modes  depending  on  the  Weber  number  as 
shown  in  Table  1. 


Table  1  Summary  of  Fuel  Spray  Modeling 


Physical  Process 

Model 

Results 

Atomizer  Free- Surface  Flow 

Ibrahim  and  Jog’s  Numerical 
Model  [13] 

Liquid  sheet  thickness 

Spray  angle 

Sheet  velocity 

Primary  Atomization 

Senecal  et  al.’s  Linearized 
Sheet  Instability  Atomization 
Model  [14] 

Ligament  size 

Droplet  size 

Secondary  Atomization 

-  Vibration  mode:  0<We<l  1 

-  Bag  mode:  1  l<We<35 

Multimode:  35<We<80 

-  Sheet- thinning:  80<We<350 

-  Catastrophic  mode:  We>350 

Hybrid  TAB/KH/RT  Model 

No  breakup:  1  l<We<3 5 

TAB  Model:  ll<We<35 
TAB/KH  Model:  35<We<80 

-  KH  Model:  80<We<350 

-  KH/RT  Model:  We>3 50 

Droplet  size 

Figure  1  Modeling  Schematics  of  High  Pressure  Fuel  Spray 


1 .  Primary  A tomization 

The  primary  atomization  of  hollow  cone  spray  is  accounted  for  by  the  LISA  model,  which  determines  the  ligament 
size  based  on  the  most  unstable  wave  length  based  on  a  linear  instability  analysis.  The  LISA  model  defines  a 
dispersion  relation  between  the  disturbance  growth  rate,  co ,  and  wave  number,  k ,  which  are  obtained  by  linearzing 
the  liquid  continuity  and  momentum  equations  as: 

co2  [tanh(^)  +  <g]  +  ft)[4i/1A:2  tanh(kh)  +  2iQkU~^  +  4vfk4  tanh(&/z) 

rrk3  (14) 

-4 v\k3  tanh  (kh)  -  QU2k 2  +  —  =  0 

A 


where  h ,  Q,  ph  U  and  are  the  liquid  sheet  thickenss,  gas  to  liquid  density  ratio,  liquid  density,  sheet  velocity  and 
liquid  viscosity  respectively.  The  dispersion  equation  is  solved  with  a  specified  liquid  sheet  thickness,  velocity  and 
fuel  properties  by  using  the  so-called  golden-search  method  to  obtain  the  maximum  instability  growth  rate  and  the 
corresponding  wavelength.  This  wavelength  is  used  to  define  the  ligament  size  as: 


dL 


(15) 


Based  on  Weber’s  column  theory,  the  droplet  size  that  breaks  away  from  the  ligament  is  determined  as: 

dD  =  l.SSdL  (1  +  30A)16 


(16) 
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Physically,  complex  multiphase  processes  are  involved  during  this  primary  atomization  event  including  the  liquid 
sheet  breakup  into  ligaments  and  ligament  breakup  into  droplets.  All  these  processes  are  computationally  simplified 
and  coupling  with  the  surrounding  gas  is  neglected  to  prevent  any  inaccurate  physical  interactions.  The  primary 
breakup  time  is  determined  empirically  as: 


x  — 


(17) 


2.  Secondary  Atomization 

Secondary  atomization  is  modeled  by  means  of  a  hybrid  model  that  combines  the  Taylor  Analogy  Breakup 
(TAB)  [15],  Reitz’s  Kelvin-Helmholtz  (KH)  and  Rayleigh-Taylor  (RT)  models  [16-18].  Depending  on  the  Weber 
number,  there  are  distinct  breakup  modes  during  the  secondary  atomization.  The  hybrid  model  adaptively  applies  the 
appropriate  submodel  corresponding  to  the  physical  breakup  mode.  Among  the  submodels,  the  TAB  model 
describes  the  bag  breakup  mode,  which  provides  good  accuracy  for  Weber  numbers  in  the  range  of  1  l<We<35.  On 
the  other  hand,  the  KH  model  represents  the  sheet-thinning  mode  (or  stripping  mode)  that  is  suitable  for 
80<We<350  and  produces  a  reasonable  number  of  broken  child  parcels,  thereby  providing  a  continuous  vaporized 
fuel.  The  transition  between  the  bag-breakup  and  sheet-thinning  modes  occurs  at  35<We<80  and  we  empirically 
determined  the  critical  Weber  number,  Wec,  to  be  40  in  the  present  simulations.  Lastly,  for  even  higher  Weber 
numbers,  RT  instability  is  found  to  occur,  which  is  driven  by  drop  acceleration  due  to  aerodynamic  drag. 


gas 


(a)  TAB  model  for  a  bag  breakup  mode:  1  l<We<35 


0 


•  *>  c 
©  •  • 

•  •  •  # 

► 

©  © 


A 


KH 


(b)  KH  model  for  a  sheet-thinning  breakup  mode:  80<We<350 


(c)  KH/RT  model  for  a  catastrophic  mode:  We>350 
Figure  2  Secondary  Atomization  Processes 


The  TAB  model  proposed  by  O’Rourke  and  Amsden  [15]  is  based  on  the  classical  analogy  between  an 
oscillating  and  distorting  droplet  and  a  spring  mass  system.  It  is  regarded  that  the  external,  restoring  and  damping 
forces  in  a  spring  mass  system  correspond  to  the  drag,  surface  tension  and  viscous  forces  in  a  drop  oscillation  and 
distortion  process.  The  TAB  model  defines  a  nondimensionalized  droplet  oscillation  parameter  to  be  y  =  x/(Cbr) 
where  x  is  the  displacement  of  the  droplet  equator  from  its  spherical  position  and  Cb  is  a  constant  equal  to  0.5.  The 
droplet  oscillation,  y,  can  be  expressed  as  Eq.  (18)  in  terms  of  fluid  properties  and  model  constants  solving  Newton’s 
2nd  law: 


y(t)  =  Wec+e~(,,‘d) 


(y0-Wec)cos(cot)  + 


CO 


dt 


o  +  JV 


-  We, 


sin(w^) 
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(18) 


The  breakup  occurs  when  y>l  and  the  Sauter  Mean  Diameter  (SMD)  after  breakup  is  determined  by  the  energy 
balance  as: 

d 


^32  ,c 


{  !  8 Ky2  |  P,r3(dy/dtf  (6K- 5 
+  20  +  cr  i  120 


(19) 

The  actual  parcel  size  in  the  simulation  is  defined  using  a  Rosin-Rommler  distribution.  Additionally  the  child  drop 
after  breakup  is  prescribed  a  velocity  component  normal  to  the  parent  drop  velocity: 


v  ,  =  C  C.r  — 
dt 


(20) 


Reitz’s  breakup  model  is  based  on  the  Kelvin-Helmholtz  instability  phenomenon  [16-18]  and  the  results  of  the 
linear  stability  analysis  are  directly  used  in  determining  the  important  parameters  in  the  atomization  process.  For 
easier  use  of  the  linear  stability  analysis  results,  the  maximum  growth  rate,  Q  ,  and  corresponding  wavelength,  A , 
of  the  Kelvin-Helmholtz  instability  mode  are  expressed  by  curve-fits  in  terms  of  nondimensional  numbers  as: 


P,a 3 

05  (0.34  +  0.38Weg2) 

G 

~  (l  +  Oh)(l  +  1.4Ta06) 

A™  902(l+°-45Oh05)(l+0-4Ta07) 


a 


(l  +  0.87Weg1S7)°6 


(21) 

(22) 


p  U-V  d  [7T 

where  Weg  is  the  Weber  number  for  the  gas  phase,  We  = - ,  Oh  is  the  Ohnesorge  number,  Oh  =  vlx  — - 

g  V  GCl 


rc  —  AAkH 

where  B0  is  a  model  constant  and  set  equal  to  0.61.  Simultaneously,  the  parent  drop  size  is  reduced  by 

da  _  ( a-r ) 

dt  t 

where  r  is  a  time  constant  and  determined  from 

r  =  3.126Bxa  /  AkhQkh 


and  Ta  is  the  Taylor  parameter,  Ta  =  Oh^Weg  .  The  liquid  breakup  is  modeled  by  adding  new  child  parcels  and 
their  size  is  determined  by 

(23) 

(24) 

(25) 

Child  parcels  are  released  when  the  stripped  mass  removed  from  the  parent  parcel  exceeds  a  few  percent  of  the 
average  injected  parcel  mass.  Newly  formed  parcels  are  assigned  a  random  velocity  direction  within  a  confined  cone 
angle  defined  by 

tan {0  /  2)  =  v41AkhQkh  /  U  (26) 

where  ^4;  =  0.188.  In  addition,  as  the  parent  parcel  reduces  in  size,  its  mass  is  preserved  by  controlling  the  drop 
numbers  contained  in  the  parcel  as  Na3  =  N0a30 . 

Lastly,  the  Rayleigh-Taylor  (RT)  model  [16]  is  based  on  aerodynamic  forces  normal  to  the  interface  between 
two  fluids  that  can  drive  the  instability.  Taking  the  external  force  acting  on  the  droplet  to  be  given  by  the 
aerodynamic  drag,  the  acceleration  can  be  expressed  by  dividing  the  drag  force  by  the  droplet  mass: 

"  pA 


a  =-CD 
8  D 


Pir 


(27) 


Based  on  the  linearized  instability  analysis  by  Chang,  the  wavelength  and  frequency  of  the  fastest  growing  waves 
are  determined  as: 


Art  =  2n  A 


(28) 


Q  =  — 

RT  V  3 


aPl 
3  G 


(29) 
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The  droplet  breakup  is  allowed  to  occur  when  the  RT  wavelength,  ART  ,  is  smaller  than  the  droplet  diameter  until  the 
RT  breakup  time,  rRT  =  1  /  QRT  . 


F.  Eulerian-Lagrangian  Coupling 

Eulerian-Lagrangian  coupling  is  accomplished  through  the  source  term  vector,  S,  in  Eq.  (2).  Mass,  chemical 
species,  momentum  and  energy  transfer  from  the  Lagrangian  phase  are  provided  as  sources  for  the  Eulerian  phase: 


dt 


r  p  ' 

/  X-1  A 

LNnmv 

pY 

Zv""1.' 

pY 

= 

Za^Ku-w/u) 

pti'-p 

-fb.Qd) 

PK  J 

{  0  J 

(30) 


where  subscript  N  and  n  represent  the  number  of  parcels  and  the  number  of  drops  in  an  individual  parcel 
respectively. 


G.  Computational  Geometry,  Mesh  and  Boundary  Conditions 

Three  cases  are  considered  in  the  present  study  at  a  given  inlet  air  temperature  and  equivalence  ratio  as 
summarized  in  Table  2.  In  the  three  simulations,  two  different  geometries  are  employed  as  shown  in  Figure  3.  The 
full  geometry  contains  all  the  major  components  in  the  combustor  as  indicated  in  Figure  3(a),  where  the  y-axis  scale 
is  exaggerated  to  clearly  illuminate  the  major  components.  The  actual  laboratory  combustor  configuration  has  a  long 
length  with  a  high  aspect  ratio  as  shown  in  Figure  3(b).  The  second  geometry,  given  in  Figure  3(c),  is  a  reduced 
geometry,  in  which  only  a  part  of  the  air  plenum  and  the  combustor  including  the  EDI  element  are  considered.  The 
main  difference  between  the  full  and  reduced  geometries  lies  in  the  ability  to  support  longitudinal  acoustic  modes  in 
the  combustor.  The  full  geometry  enables  this  by  isolating  the  acoustic  waves  in  the  combustor  by  the  insertion  of  a 
slotted  inlet  and  an  exit  nozzle,  whereas  the  reduced  geometry  is  acoustically  open. 


Table  2  Summary  of  Cases  Considered 


Fabel 

Tair 

0 

Geometry 

Chemistry 

Case  1 

700  K 

0.60 

Reduced/ Acoustically  Open  (Fig.  3  (c)) 

No  Spray  and  Reaction 

Case  2 

Reduced/ Acoustically  Open  (Fig.  3  (c)) 

Spray  Combustion 

Case  3 

Full/ Acoustically  Closed  (Fig.  3  (b)) 

Spray  Combustion 

In  the  simulation,  all  the  surfaces  of  the  geometry  except  for  the  inlet  and  exit  planes  are  treated  as  no-slip  and 
adiabatic  walls.  The  inlet  plane  is  located  at  the  top  left  end  of  Figure  3  for  the  full  geometry  and  left-end  surface  for 
the  reduced  geometry.  The  air  flow  enters  through  the  inlet  using  a  constant  mass  flow  condition.  At  the  exit  nozzle 
located  at  the  right-end  surface,  an  outlet  condition  is  imposed  by  a  characteristic  back  pressure  condition.  The  fuel 
spray  is  injected  at  the  atomizer  exit  plane  and  the  number  of  drops  is  specified  in  order  to  account  for  the  correct 
fuel  mass  flow  rate. 

A  hexahedral  structured  mesh  system  is  used  for  the  geometry  as  shown  in  Figure  4  and  “butterfly”  meshing  is 
used  for  the  circular  cross-sections.  A  total  of  six  million  cells  are  employed  for  the  full  geometry  and  one  million 
cells  are  used  for  the  reduced  geometry.  According  to  our  grid  refinement  studies  using  the  full  geometry,  coarser 
meshes  of  one  to  two  million  cells  are  found  to  be  incapable  of  properly  resolving  the  acoustic  oscillations. 
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Combustor 


4- 


L, 


^4 


U 


(a)  Full  combustor  highlighting  major  components 

I"  c  > 

(b)  Full  Combustor  Geometry  with  1 :2  aspect  ratio 


(c)  Acoustically-Open  Combustor  Geometry  with  1 :2  aspect  ratio 

Figure  3  Computational  Domain  of  the  laboratory  LDI  model  combustor 


Figure  4  Computational  Mesh  visualized  around  the  fuel  nozzle  and  swirler 


III.  Decomposition  Methods  For  Combustion  Dynamics  Diagnostics 

To  understand  the  underlying  structure  of  the  acoustic  modes  and  the  associated  combustion  dynamics,  we  utilize 
Dynamic  Mode  Decomposition  (DMD)  techniques.  DMD  takes  an  ensemble  of  data,  either  instantaneous  or  a  series 
of  snapshots,  as  inputs  and  uses  the  Arnoldi  algorithm  (for  DMD)  to  define  ‘modes’  that  are  mathematically 
orthogonal  to  each  other.  A  detailed  mathematical  description  of  the  DMD  approach  can  be  found  in  Rowley  et  al. 
[19]  and  Schmid  [20].  Here,  the  mathematical  model  is  briefly  reviewed. 

A.  Mathematical  Model  of  DMD 

We  wish  to  approximate  a  function  z(x,t)  over  some  domain  of  interest  as  a  finite  sum  in  the  variables-separated 
form 

M 

(31) 

k= 1 
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with  a  reasonable  expectation  that  the  approximation  becomes  exact  as  M  approaches  infinity,  except  possibly  on  a 
set  of  measure  zero.  Note  that  in  Eq.  (31)  there  is  no  fundamental  difference  between  x  and  t ,  but  we  usually  think 
of  X  as  a  spatial  coordinate  (ID,  2D  &  3D  vector- valued)  and  t  as  the  temporal  coordinate. 

The  representation  of  Eq.  (31)  is  not  unique.  For  example,  if  the  domain  (experimental  or  computational)  is  a 
bounded  interval  X  on  the  real  line,  then  the  functions  0^(x)can  be  chosen  as  a  Fourier  series,  or  Legendre 

polynomials,  or  Chebyshev  polynomials,  and  so  on.  For  different  selections  of  the  space-dependent  function,  0*(x) 

,  the  corresponding  time-dependent  function,  ak(t)  ,  will  be  different.  They  can  be  periodic  or  non-periodic,  single¬ 
frequency  dominated  or  multi-frequency  dominated. 

In  the  POD  analysis,  the  spatial  functions,  0*(x) ,  are  chosen  to  be  orthonormal  functions,  i.e. 


j  (*)$*,  (x)dx 

X 


1  if  kx  =k2 
0  otherwise 


(32) 


Then 


ak(t)  =  J  z(x,ty&k(x)dx,  (33) 

x 

From  Eq.  (32)  and  (33),  it  is  noted  that  given  that  0*(x)  is  selected  to  be  orthonormal,  the  determination  of  ak(t) 
will  be  only  dependent  on  0*(x)  rather  than  on  the  other  O  s. 


B.  Arnoldi  Algorithm 

In  DMD  analysis,  in  order  to  obtain  single  frequency  dynamic  modes,  linear  mapping  is  assumed  from  one 
snapshot  to  another.  Suppose  the  dataset  is  represented  as  a  snapshot  sequence, 


V{  ={vpv2,v3, . ,v*} 


(34) 


where  v,.  stands  for  the  ilh  flow  field  and  vj+1  =  Bvi .  The  matrix  B  here  represents  the  linear  mapping  matrix. 
Therefore, 


V1N={  v19Bv19B2v19 9Bn~\} 


(35) 


Another  assumption  is  that  there  exists  a  specific  number  A,  beyond  which  the  vector  can  be  expressed  as  linear 
combination  of  the  previous  vectors, 

Vj  =  axvx  +a2v2  +  ...  +  aN_xvN_x  or  vx  =  VxN~la  +  r  (36) 

Hence, 

BV*-1  =  V2  =  VxN~lS  +  reTN_x  (37) 

where. 


S  = 


0 

1  0 


10  a 


N-2 


1  a 


N- 1  J 


Applying  the  eigen-value  decomposition  for  matrix  S, 


S  =  T~lAT 


(38) 


(39) 


where  matrix  T  is  the  eigen-vector  matrix  of  S.  Suppose  a  sufficient  number  of  snapshots  are  used,  the  eigenvalues 
of  S  can  very  well  represent  the  eigenvalues  of  A,  which  contains  the  time-evolution  information  of  the  flow  field. 
Also,  the  dynamic  modes  corresponding  to  single  frequency  response  can  be  constructed  as, 

O  =  VX~XT~X  (40) 

Similarly,  the  original  data  set  can  be  decomposed  into  the  form  in  Eq.  (31), 


vx  =  or 


(41) 


where  matrix  O  contains  the  dynamic  spatial  information  and  matrix  T  contains  the  temporal  evolutional 
information. 
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IV.  Results 

The  diverse  physics  in  the  combustor  and  their  coupling  mechanisms  are  investigated  by  a  series  of  combustion 
dynamics  simulations,  which  illuminate  the  roles  of  hydrodynamics,  spray  formation,  heat  release,  and  acoustics. 
The  cases  listed  in  Table  2  are  designed  to  focus  on  the  specific  flow  physics  of  interest.  For  example,  Cases  1  and  2 
show  how  the  flow  characteristics  are  influenced  by  chemical  reactions,  especially  the  combustion  heat  release.  On 
the  other  hand,  Case  3  isolates  the  acoustics  effects  that  are  the  result  of  the  acoustically  closed  environment.  In  all 
cases,  the  acoustic  characteristics  from  the  computational  results  are  compared  with  those  from  the  companion 
experimental  study. 

A.  Dynamic  Flow  Characteristics 

A  large  central  toroidal  recirculation  zone  (CTRZ)  is  found  in  time-averaged  results  in  both  the  non-reacting  and 
reacting  flow  cases  as  shown  in  Figure  5.  The  CTRZ  shows  the  presence  of  asymmetric,  unsteady  spiral  motion. 
This  behavior  is  commonly  found  in  swirling  flows  and  is  associated  with  self-excited  hydrodynamic  flow 
instabilities  referred  to  as  precessing  vortex  core  (PVC)  instabilities.  The  strong  swirling  flow  generated  by  the 
swirler  vane  has  radial  gradients  of  the  azimuthal  velocity  and  pressure.  As  the  flow  is  expanded  continuously  in  the 
venturi-diverging  section  and  subsequently  at  the  combustor  backstep,  there  is  a  corresponding  decay  of  the 
azimuthal  velocity  and  the  radial  pressure  gradients.  In  turn,  this  leads  to  a  reversal  in  the  axial  pressure  gradient 
which  results  in  the  formation  of  the  CTRZ.  Since  the  axial  decay  of  azimuthal  momentum  and  radial  pressure 
gradient  forms  the  CTRZ  structure,  the  combustor  head  geometry  and  flow  conditions  can  be  critical  factors  in 
determining  the  size  and  strength  of  CTRZ.  The  PVC  instabilities  occur  when  the  CTRZ  becomes  unstable  and 
asymmetric.  A  series  of  self-excited  vortices  is  apparent  in  the  planar  view  of  the  combustor  head  region  shown  in 
Figure  6.  The  vortices  are  generated  at  the  nozzle  exit  and  merge  into  the  main  recirculation  zone  (frames  b  to  d). 
This  process  alternates  during  the  PVC  time  period,  inducing  an  asymmetric  vortex  structure. 


Figure  5  Time-averaged  streamlines:  non-reacting  flow  (left)  and  reacting  flow  (right) 
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(a)  t  =  t0 


(b)  t  —  to  +  1/6  ipvc  (c)  t  —  to  +  2/6  ipvc 


(d)  t  —  to  +  3/6  ipvc 


(e)  t  —  to  +  4/6  ipvc 


(f)  t  —  to  +  5/6  ipvc 


Figure  6  Instantaneous  streamlines  during  a  PVC  mode  period 

The  main  difference  between  non-reacting  and  reacting  flows  is  the  existence  of  a  secondary  downstream 
vortical  structure  in  the  non-reacting  case.  The  volume  expansion  of  the  gas  due  to  combustion  allows  the  main 
recirculation  zone  to  fill  a  larger  region  at  the  chamber  head.  The  main  recirculation  zone  then  serves  as  a  blockage 
for  the  swirling  air  flow,  the  mean  flow  speed  at  the  outer  shear  layer  is  accelerated.  This  increased  axial  momentum 
in  turn  leads  to  the  straight  flow-through  pattern  observed  beyond  the  CTRZ  in  the  reacting  flow  case.  In  addition, 
the  expansion  of  CTRZ  and  the  acceleration  of  the  air  flow  also  leads  to  a  shifting  of  the  PVC  frequency  to  higher 
values.  The  distinct  frequencies  are  obtained  using  the  DMD  method  and  are  summarized  in  Table  4,  indicating  a 
frequency  shift  of  about  15%  between  Case  1  (non-reacting)  and  Case  2  (reacting). 

Table  3  Comparison  of  hydrodynamic  modes  for  non-reacting  and  reacting  flows 


Hydrodynamic  Mode 

Non-reacting  Flow  (Hz) 

Reacting  Flow  (Hz) 

PVC  (fundamental) 

3049 

3505 

PVC  (2nd  harmonic) 

6177 

7075 

PVC  (3rd  harmonic) 

9304 

10580 

PVC  (4th  harmonic) 

12353 

14213 

B.  Spray  and  combustion  characteristics 

In  the  present  combustion  dynamics  simulation,  the  flame  seems  to  be  strongly  influenced  by  the  hydrodynamics 
of  the  swirling  air  flow.  In  Figure  7,  which  shows  the  time-averaged  heat  release  contours,  the  flame  starts  from  the 
stagnation  point  of  CTRZ  and  follows  the  outer  shear  layer  until  the  combustor  head.  Further,  the  flame  surface  also 
fluctuates  due  to  the  PVC  instabilities.  A  series  of  shed  vortices  distorts  the  outer  shear  layer  and  the  resulting  flame. 

The  fuel  spray  is  finely  atomized  in  the  diverging  section  of  the  venturi  as  shown  in  Figure  8,  with  a  D32  ranging 
from  30  to  40  pm.  The  high  pressure  condition  rapidly  promotes  secondary  atomization  and  the  subsequent 
vaporization  rate  of  droplets.  It  is  important  to  note  that  the  strong  intensity  of  heat  release  is  enhanced  by  the 
secondary  atomization  and  is  potentially  an  important  mechanism  for  the  excitation  of  combustion  instability. 
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Figure  8  Instantaneous  spray  drops  colored  by  size 


C.  Modal  Analysis  in  an  Acoustically-Open  Combustor 

Modal  analysis  was  conducted  using  PSD’s  based  on  the  wall  pressure  at  the  chamber  head  as  shown  in  Figure 
9.  It  can  be  observed  that  the  pressure  oscillation  levels  are  amplified  in  the  presence  of  chemical  reactions.  The 
pressure  amplitude  is  0.2  kPa  for  the  non-reacting  flow  case,  whereas  it  is  0.6  kPa  for  the  reacting  flow  case.  The 
maximum  oscillation  amplitude  was  6  kPa  around  the  LDI  element  where  the  PVC  instability  occurs.  In  fact,  the 
PVC  instability  appears  to  be  a  key  driver  of  the  pressure  oscillations  in  both  cases.  In  the  non-reacting  case, 
although  the  PVC  instability  has  no  direct  mechanism  for  generating  the  acoustic  perturbations,  the  results  indicate  a 
weak  intensity  pressure  perturbations  corresponding  to  the  PVC  modes.  In  the  reacting  flow  case,  there  are  several 
distinct  pressure  peaks  in  the  PSD  plot,  which  can  be  further  investigated  using  DMD  analysis. 


(a)  non-reacting  flow 
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(b)  reacting  flow 

Figure  9  Wall  Pressure  PSDs  for  acoustically-open  combustor 


Using  the  DMD  analysis,  the  mode  power  spectrum  based  on  the  pressure  over  the  entire  flow  field  is  obtained 
as  shown  in  Figure  10.  Strong  harmonic  PVC  peaks  are  dominant  features  in  the  non-reacting  flow  case.  Additional 
pressure  peaks  appear  in  the  mode  power  spectrum  in  the  reacting  flow  case.  Longitudinal  acoustics  are  involved 
even  in  an  acoustically-open  chamber  geometry  because  the  exit  boundary  is  not  a  perfect  non-reflecting  boundary 
condition.  Moreover,  a  tangential  acoustic  mode  is  apparent  at  9433  Hz.  In  addition,  a  spinning  mode  and  associated 
higher  harmonics  are  observed  at  510,  1020  and  1530  Hz  in  the  DMD  spectrum  for  the  azimuthal  vorticity,  as  seen 
in  Figure  11.  The  major  mode  shapes  obtained  from  the  DMD  analysis  are  illustrated  in  Figure  12.  The  results 
clearly  indicate  that  the  flows  at  the  combustor  head  and  in  the  venturi  diverging  section  are  strongly  influenced  by 
the  PVC  instability  and  the  swirl  spinning  mode. 


(b)  reacting  flow 

Figure  10  DMD  Mode  Power  Spectrum  for  Pressure 


Figure  11  DMD  Mode  Power  Spectrum  for  Azimuthal  Vorticity  in  Reacting  Flow 
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Table  4  Summary  of  DMD  Modes  and  Corresponding  Frequencies 


DMD  Mode 

Frequency  (Hz) 

Unknown  mode 

892 

Spinning  mode  (1,  2,  and  3) 

510,  1020  and  1530 

PVC  mode  (1,  2  and  3) 

3505,  7075  and  10580 

Longitudinal  mode  (1L,  2L  and  3L) 

2486,  3952  and  5545 

Tangential  mode  (IT) 

9433 

c)  1L  d)lT 

Figure  12  DMD  Mode  Shapes  for  Pressure  in  an  Acoustically-Open  Chamber 


D.  Estimation  of  Self-Excited  Combustion  Instability 

Self-excited  combustion  instabilities  are  observed  in  the  acoustically-closed  combustor  flow  simulations.  In  thic 
ase,  the  present  simulations  are  preliminary  in  nature  and  did  not  include  the  full  atomization  models  as  described  in 
Section  2.E.  Instead,  the  primary  atomization  model  was  neglected  and  the  injected  drop  diameter  was  assumed  to 
be  equal  to  the  liquid  sheet  thickness  obtained  from  the  Eulerian  VOF  simulation  for  the  atomizer  free  surface  flow. 
Moreover,  the  secondary  atomization  model  was  restricted  to  the  Kelvin-Helmholtz  model.  Due  to  the  small  injected 
droplet  sizes,  secondary  atomization  effects  are  not  strong.  However,  as  indicated  in  Figure  13,  the  qualitative  trends 
in  the  PSD  plots  are  consistent  with  the  experimental  results.  Both  simulations  and  experiments  show  evidence  of  a 
strong  4L  mode  and  a  second  6000  Hz  peak  as  discussed  in  the  next  section. 


Frequency  [Hz]  Frequency  [Hz] 

(a)  Computation  (b)  Experiments 

Figure  13  Wall  Pressure  PSDs  for  acoustically-closed  combustor 
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E.  Nonlinear  PVC-Acoustic  Couplings  and  Feedback  Loop 


The  dominat  4L  and  6000  Hz  modes  observed  in  the  acoustically  closed  results  arise  from  non-linear  interactions 
between  hydrodynamics  and  acoustics  modes.  Experiments  in  lean  premixed  swirl  combustors  have  also  indicated 
the  presence  of  such  interactions.  A  hypothesis  regarding  nonlinear  PVC-acoustics  coupling  can  be  developed  using 
a  Rayleigh  criterion-based  analysis  in  the  frequency  domain.  Our  approach  uses  the  acoustic  energy  equation  written 
as: 


d 

+ 

a ,  ,  ,x 

+ — (pu  )  = 

^  pQ 

dt 

_2 pcz  2  j 

dxy  J 

_  Y 

(42) 


where,  if  we  note  that  the  pressure  perturbation  typically  originates  from  the  acoustics  and  while  the  unsteady  heat 
release  is  largely  influenced  by  the  PVC  modes.  The  pressure  and  heat  release  perturbations  can  therefore  be 
expressed  as: 

p'  -  pelo>at  □  pcos(coj)  (43) 

Q'  =  Qeia w  □  Q  cos  (a>PVCt  +  </>)  (44) 

It  is  important  to  note  the  nonlinear  pressure-heat  release  coupling  term  on  the  right  hand  side  of  Eq 
leads  to  two  interactive  frequencies  which  can  be  identified  as: 

P'Q'  =  pQcOs(a>at)cOs(G)pvct  +  </>) 

=  ^pjcos  \_(coa  +^PFC)^  +  ^]  +  COs[(ft)a  “  (°pvc)t  -  ^]} 

The  sum  and  minus  of  two  driving  frequencies  can  thus  be  produced  by  this  nonlinear  pressure-heat  release 
coupling  process.  In  the  LDI  combustor,  the  two  driving  frequencies  are  the  acoustics  and  hydrodynamics,  which  in 
turn  produce  two  interactive  modes. 

Another  noteworthy  result  is  the  establishment  of  a  feedback  loop  arising  from  this  nonlinear  interaction.  The 
interactive  modes  can  the  couple  with  other  existing  modes  such  as  the  longitudinal  and  tangential  acoustic  modes 
and  the  PVC  modes.  If  we  consider  PVC-4L  (3200  Hz/1600  Hz)  coupling,  the  interactive  modes  are  4800  Hz  and 
1600  Hz.  1600  Hz  is  the  4L  mode  and  the  PVC-4L  interaction  leads  to  further  4L  mode  amplification.  In  addition  to 
PVC-4L  coupling,  a  strong  PVC- IT  (3200  Hz  -  9600  Hz)  coupling  can  produce  interactive  modes  at  12800  Hz  and 
6400  Hz.  Here,  the  6400  Hz  mode  is  the  second  harmonic  of  PVC  which  can  also  form  a  feedback  loop,  too.  Based 
on  this  approach,  we  can  estimate  that  the  1600  Hz  (4L,  PVC-4L)  and  6400  Hz  (2PVC,  PVC- IT)  are  the  dominant 
interaction  modes,  which  agrees  well  with  strong  peaks  observed  in  the  PSD  in  Figure  13. 

In  summary,  the  combustion  dynamics  mechanism  in  the  LDI  combustor  can  be  attributed  to  nonlinear 
interactions  between  the  hydrodynamics,  chamber  acoustics,  and  combustion  as  shown  in  Figure  14.  Strong 
hydrodynamic  modes  such  as  the  PVC  instabilities  can  influence  the  flame  dynamics  and,  in  turn,  the  acoustics.  The 
unsteady  heat  release  from  the  flame  and  pressure  perturbation  from  the  acoustics  produces  two  interactive  modes. 
When  these  interactive  modes  are  associated  with  existing  hydrodynamic  or  acoustic  modes,  a  feedback  loop  is 
established  which  can  give  rise  to  the  excitation  of  combustion  instability. 


.  (42).  This  term 

(45) 


Figure  14  Schematics  of  Combustion  Dynamics  Mechanism  in  a  LDI  Combustor 
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V.  Conclusions 

Self-excited  combustion  instabilities  are  computationally  investigated  for  the  LDI  gas  turbine  combustor.  Three 
simulations  are  considered  in  the  present  investigation:  non-reacting  and  reacting  flows  in  an  acoustically  open 
chamber  and  reacting  flow  in  an  acoustically  closed  chamber.  Findings  from  simulations  can  be  summarized  as: 

•  Dynamic  flows  appear  in  both  non-reacting  and  reacting  flows.  PVC  instabilities  represent  key 
hydrodynamics  phenomena  in  the  LDI  combustor  and  the  spray  and  the  flame  respond  strongly  to  these  instabilities. 

•  Self-excited  combustion  instabilities  are  captured  successfully  in  the  simulations  and  the  dominant  peaks  in 
PSD  agree  well  with  those  in  the  experimental  results. 

•  A  hypothesis  about  the  underlying  combustion  dynamics  mechanism  in  a  LDI  gas  turbine  combustor  has 
been  presented,  wherein  the  key  driver  is  the  nonlinear  interaction  between  the  P VC-related  heat-release  modes  and 
the  acoustics  modes.  We  argue  that  the  presence  of  a  feedback  loop  may  lead  to  the  excitation  and  sustenance  of 
combustion  instabilities. 
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